## RTMB dat and par setup ####
# Dat
dat <- list(srdat = srdat,
WAbase = WAbase,
WAin = WAin,
lineWA = seq(min(WAbase$logWAshifted),
max(WAbase$logWAshifted), 0.1), # Not added to NLL
mean_logWA = mean_logWA,
logRS = log(srdat$Rec) - log(srdat$Sp),
prioronly = 0) # 0-run with data, 1-prior prediction mode
# External vectors
N_Stk <- max(srdat$Stocknumber + 1)
N_Obs <- nrow(srdat)
stk = srdat$Stocknumber + 1
lhdiston <- T
bias.cor <- F
# Parameters/Initial values
par <- list(b0 = c(10, 0), # Initial values for WA regression intercepts
bWA = c(0, 0), # Inital values for WA regression slopes
logSREP_re = numeric(N_Stk), # Zeroes
logAlpha0 = 0.6,
logAlpha_re = numeric(nrow(dat$WAbase)), # Zeroes
tauobs = 0.01 + numeric(N_Stk), # Constrained positive
logSREP_sd = 1,
logAlpha_sd = 1
)
if (lhdiston) {
par$logAlpha02 <- 0
}
f_srep <- function(par){
getAll(dat, par)
N_Stk = max(srdat$Stocknumber + 1) # number of stocks
stk = srdat$Stocknumber + 1 # vector of stocknumbers
N_Obs = nrow(srdat) # number of observations
N_Pred = nrow(WAin) # number of predicted watershed areas
S = srdat$Sp
type = lifehist$lh
type_tar = as.numeric(WAin$lh)
SREP <- numeric(N_Stk)
# logE_pred <- numeric(N_Stk)
logSREP <- numeric(N_Stk)
logAlpha <- numeric(N_Stk)
logRS_pred <- numeric(N_Obs)
SREP_tar <- numeric(N_Pred)
logSREP_tar <- numeric(N_Pred)
logAlpha_tar <- numeric(N_Pred)
# Simulated line vectors
line <- length(lineWA)
logSREP_line_stream <- numeric(line)
SREP_line_stream <- numeric(line)
logSREP_line_ocean <- numeric(line)
SREP_line_ocean <- numeric(line)
if (bias.cor) {
biaslogSREP <- -0.5*logSREP_sd^2
biaslogAlpha <- -0.5*logAlpha_sd^2
biaslogRS <- -0.5*(sqrt(1/tauobs))^2
} else {
biaslogSREP <- 0
biaslogAlpha <- 0
biaslogRS <- numeric(N_Stk)
}
nll <- 0 # Begin negative log-likelihood
nll <- nll - sum(dnorm(b0[1], 10, sd = 31.6, log = TRUE)) # Prior
nll <- nll - sum(dnorm(b0[2], 0, sd = 31.6, log = TRUE)) # Prior
nll <- nll - sum(dnorm(bWA[1], 0, sd = 31.6, log = TRUE)) # Prior
nll <- nll - sum(dnorm(bWA[2], 0, sd = 31.6, log = TRUE)) # Prior
nll <- nll - sum(dnorm(logAlpha0, 0.6, sd = 0.45, log = TRUE)) # Prior (rM)
if(lhdiston) nll <- nll - sum(dnorm(logAlpha02, 0, sd = 31.6, log = TRUE)) # Prior (rD)
## Second level of hierarchy - Ricker parameters:
for (i in 1:N_Stk){
nll <- nll - dnorm(logSREP_re[i], 0, sd = 1, log = TRUE)
logSREP[i] <- b0[1] + b0[2]*type[i] + (bWA[1] + bWA[2]*type[i]) * WAbase$logWAshifted[i] + logSREP_re[i]*logSREP_sd + biaslogSREP
SREP[i] <- exp(logSREP[i])
nll <- nll - dnorm(logAlpha_re[i], 0, sd = 1, log = TRUE)
if(lhdiston) logAlpha[i] <- logAlpha0 + logAlpha02*type[i] + logAlpha_re[i]*logAlpha_sd + biaslogAlpha
else logAlpha[i] <- logAlpha0 + logAlpha_re[i]*logAlpha_sd + biaslogAlpha
nll <- nll - dgamma(tauobs[i], shape = 0.0001, scale = 1/0.0001, log = TRUE)
}
## First level of hierarchy: Ricker model:
for (i in 1:N_Obs){
logRS_pred[i] <- logAlpha[stk[i]]*(1 - S[i]/SREP[stk[i]]) + biaslogRS[stk[i]]
if(!prioronly){
nll <- nll - dnorm(logRS[i], logRS_pred[i], sd = sqrt(1/tauobs[stk[i]]), log = TRUE)
}
}
## Calculate SMSY for Synoptic set - for plotting
SMSY_r = numeric(nrow(WAbase))
BETA_r = numeric(nrow(WAbase))
for (i in 1:N_Stk){
BETA_r[i] <- logAlpha[i] / SREP[i]
SMSY_r[i] <- (1 - LambertW0(exp(1 - logAlpha[i]))) / BETA_r[i]
}
## PREDICTIONS
BETA = numeric(nrow(WAin))
SMSY = numeric(nrow(WAin))
SGEN = numeric(nrow(WAin))
for (i in 1:N_Pred){
if(lhdiston) logAlpha_tar[i] <- logAlpha0 + logAlpha02*type_tar[i] + biaslogAlpha
else logAlpha_tar[i] <- logAlpha0 + biaslogAlpha
logSREP_tar[i] <- b0[1] + b0[2]*type_tar[i] + (bWA[1] + bWA[2]*type_tar[i])*WAin$logWAshifted_t[i] + biaslogSREP
SREP_tar[i] <- exp(logSREP_tar[i])
# Predict BETA
BETA[i] <- logAlpha_tar[i]/SREP_tar[i]
# Predict SMSY
SMSY[i] <- (1-LambertW0(exp(1-logAlpha_tar[i])))/BETA[i]
# Predict SGEN
SGEN[i] <- -1/BETA[i]*LambertW0(-BETA[i]*SMSY[i]/(exp(logAlpha_tar[i])))
}
# Create predictions on an simulated line
for (i in 1:line){
logSREP_line_ocean[i] <- b0[1] + b0[2] + (bWA[1] + bWA[2])*lineWA[i] + biaslogSREP # DATA - not in likelihood
SREP_line_ocean[i] <- exp(logSREP_line_ocean[i])
logSREP_line_stream[i] <- b0[1] + (bWA[1])*lineWA[i] + biaslogSREP # DATA - not in likelihood
SREP_line_stream[i] <- exp(logSREP_line_stream[i])
}
## ADREPORT - internal values (synoptic specific/Ricker)
REPORT(b0) # Testing simulate()
REPORT(bWA) # Testing simulate()
REPORT(logRS_pred)
alpha <- exp(logAlpha)
# REPORT(logRS) # logRS for all 501 data points
REPORT(logSREP_re)
REPORT(logSREP_sd)
REPORT(SREP) # E (Srep) for all synoptic data set rivers (25)
REPORT(logSREP)
REPORT(logAlpha) # model logAlpha (25)
REPORT(logAlpha0)
REPORT(logAlpha02)
REPORT(logAlpha_re) # random effect parameter for resampling
REPORT(logAlpha_sd)
REPORT(alpha)
REPORT(SMSY_r)
REPORT(BETA_r)
REPORT(tauobs) # Necessary to add back in observation error?
alpha_tar <- exp(logAlpha_tar)
REPORT(SREP_tar)
REPORT(logSREP_tar)
REPORT(logAlpha_tar)
REPORT(alpha_tar)
REPORT(BETA)
REPORT(SMSY)
REPORT(SGEN)
# Simulated line values for plotting
REPORT(SREP_line_stream)
REPORT(logSREP_line_stream)
REPORT(SREP_line_ocean)
REPORT(logSREP_line_ocean)
nll # output of negative log-likelihood
}